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Abstract We study a conservative 5-point cell-centered 
finite volume discretization of the high-contrast diffu- 
sion equation. We aim to construct preconditioners that 
are robust with respect to the magnitude of the coeffi- 
cient contrast and the mesh size simultaneously. For 
that, we prove and numerically demonstrate the ro- 
bustness of the preconditioner proposed by Aksoylu et 
al. (2008, Comput. Vis. Sci. 11, pp. 319-331) by ex- 
tending the devised singular perturbation analysis from 
linear finite element discretization to the above dis- 
cretization. The singular perturbation analysis is more 
involved than that of finite element because all the sub- 
blocks in the discretization matrix depend on the diffu- 
sion coefficient. However, that dependence is eliminated 
asymptotically. This allows the same preconditioner to 
be utilized due to similar limiting behaviours of the sub- 
matrices; leading to a narrowing family of precondition- 
ers that can be used for different discretizations — a de- 
sirable preconditioner design goal. We compare our nu- 
merical results to standard cell-centered multigrid and 
observe that performance of our preconditioner is in- 
dependent of the utilized prolongation operators and 
smoothers. 

As a side result, we also prove that the solution over 
the highly-diffusive island becomes constant asymptot- 
ically. Integration of this qualitative understanding of 
the underlying PDE to our preconditioner is the main 
reason behind its superior performance. Diagonal scal- 
ing is probably the most basic preconditioner for high- 
contrast coefficients. Extending the matrix entry based 
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spectral analysis introduced by Graham and Hagger, we 
rigorously show that the number of small eigenvalues of 
the diagonally scaled matrix depends on the number of 
isolated islands comprising the highly-diffusive region. 
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1 Introduction 

We advocate that qualitative understanding of the PDE 
operators and their dependence on the coefficients is 
essential for designing preconditioners. Since, the per- 
formance, hence the robustness, of a preconditioner de- 
pends essentially on the degree to which the precon- 
ditioner approximates the properties of the underlying 
PDE. Therefore, designing preconditioners involves a 
process that draws heavily upon effective utilization of 
tools from operator theory as well as singular pertur- 
bation analysis (SPA). In the operator theory frame- 
work, Aksoylu and Beyer |31[S] have studied the diffu- 
sion equation with rough coefficients. The roughness of 
coefficients creates serious complications. For instance, 
it was shown in [4 that the standard elliptic regularity 
in the smooth coefficient case fails to hold. Moreover, 
the domain of the diffusion operator heavily depends 
on the regularity of the coefficients. 

Roughness of PDE coefficients causes loss of robust- 
ness of preconditioners. We aim to establish robustness 
with respect to the magnitude of the coefficient con- 
trast and the mesh size simultaneously. In that regard, 
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SPA provides valuable insight into the qualitative na- 
ture of the underlying PDE. In the case of linear finite 
element (FE), Aksoylu et al. [7] devised a SPA on the 
matrix entries to study the robustness of the same pre- 
conditioner Backs under consideration in this article. 
SPA turned out to be an effective tool in analyzing 
certain behaviors of the discretization matrix K such 
as the asymptotic rank, decoupling, low-rank perturba- 
tions of the resulting submatrices. This information in 
turn is exploited to accomplish dramatic computational 
savings. In [7], we also provided a rigorous convergence 
analysis of Bagks- Hence, with the insights provided 
by SPA and the operator theory, we are in control of the 
effectiveness and computational feasibility at the same 
time. 

The preconditioner Backs originates from the fam- 
ily of robust preconditioners constructed by Aksoylu 
and Klie [5] for the cell-centered finite volume (FV) dis- 
cretization of the high-contrast diffusion equation for 
porous media flow related applications based on the 
two-point approximation scheme studied in [2]. How- 
ever, for the original variants of BackSi robustness 
with respect to the contrast size was the main design 
feature. Hence, the emphasis in [H] was mostly on de- 
flation methods (that are used as stage-two precondi- 
tioners) based on the Krylov subspace solvers. As in 7 , 
in this article, we incorporate (cell-centered) multigrid 
preconditioners to restore robustness with respect to 
the mesh size. Hence, the main purpose of this arti- 
cle is to extend the preconditioner Backs f^nd the re- 
lated analysis to a 5-point conservative cell-centered FV 
discretization. From the flow application perspective, 
maintaining the continuity of the flux across the con- 
trol volume interfaces ensures the highly popular and 
crucial discretization property of local mass conserva- 
tion. Our interest in flow applications is the main rea- 
son behind conducting research in the direction of mass 
conservative discretizations. 

We prove and numerically demonstrate that the very 
same preconditioner Backs that is used for FE dis- 
cretization can also be used with minimal modification 
for FV discretizations. This was possible by the help of 
SPA because we have identified similar algebraic fea- 
tures of the discretization matrices between linear FE 
and FV methods. The same preconditioner can be uti- 
lized due to similar limiting behaviours of the subma- 
trices. This observation leads to a narrowing family of 
preconditioners that can be used for different discretiza- 
tions. Therefore, we have accomplished to construct a 
preconditioner that is compatible with and equally ef- 
fective under different discretizations, and this is a de- 
sirable feature in the design and construction of pre- 
conditioners. In addition, extension to FV discretiza- 



tion does not spoil Backs^^ inherit algebraic nature. 
The first algebraic phase involves partitioning of the 
degrees of freedom (DOF) into a set corresponding to 
a high-diffusivity and a low-diffusivity region. For high 
enough contrast, we can still obtain the partitioning by 
examining the diagonal entries of K. This simple al- 
gebraic examination rules out the standard multigrid 
requirement of the coarsest mesh resolving the bound- 
ary of the island. 

The diffusion equation with discontinuous coefficients 
is known as the interface problem in the computational 
ffuid dynamics community. There has been intense re- 
search activity on the interface problem. It is such a 
well-established problem that one can find it in the text 
books dedicated to multigrid methods, for instance, by 
Wesseling ^\ and Trottenberg, Oosterlee, and Schiiller 
[25] . Mohr and Wienands [21 revisited the cell-centered 
multigrid (CCMG) preconditioner and attributed the 
pioneering CCMG to Wesseling [29]; see [21] and the 
references therein for further review of CCMG. 

For interface problems, there have been many at- 
tempts to construct problem independent prolongation 
and restriction operators to accommodate the the rough- 
ness of coefficients. Among the variational approaches, 
Wesseling and Wesseling and Khalil constructed 
such prolongations for high-contrast coefficients, whereas 
Kwak [T7] studied medium-contrasts. Kwak and Lee [TS] 
proposed problem dependent prolongations for medium- 
contrasts. Among the non-variational approaches, Ew- 
ing and Shen [9] examined medium-contrast with piece- 
wise constant prolongation and restriction operators. 
In our CCMG and Backs implementations, we prefer 
problem independent prolongation operators utilizing 
Galerkin variational approach based on Wesseling and 
Khalil [3T] and bi- linear interpolation; see ^1] Figure 
2] or P- 72]. 

Diagonal scaling is probably the most basic pre- 
conditioner for discretizations with high-contrast coef- 
ficients. Although diagonal scaling has no effect on the 
asymptotic behaviour of the condition number, Gra- 
ham and Hagger [TH] observed, in the case of FE, that 
spectrum of the diagonally scaled matrix A enjoys bet- 
ter clustering than that of K. The spectrum of A is 
bounded from above and below except a single eigen- 
value in the case of a single isolated highly-diffusive 
island. On the other hand, the spectrum of K con- 
tains eigenvalues approaching infinity with cardinality 
depending on the number of DOF contained within 
highly-diffusive island. For faster convergence, the Krylov 
subspace solver favors the clustering provided by diago- 
nal scaling. Based the matrix entry based spectral anal- 
ysis introduced by Graham and Hagger [T3] for linear 
FE, we extend the spectral analysis to cell-centered FV 
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discretization and obtain a similar spectral result for A. 
Namely, the number of small eigenvalues of A depends 
on the number of isolated islands comprising the highly- 
diffusive region. This rigorous analysis is presented in 
Section g] 

We extend the devised SPA from FE to FV dis- 
cretization in order to explain the properties of the 
submatrices related to K{m). In particular, SPA of 
KHHirn), as diffusivity m ^ oo, has important im- 
plications for the behaviour of the Schur complement 
S{ni) of KffHi'nT') in K{m). Namely, 

S{m) = 5*°° + ©(m^i) , 

where S°° is a low rank perturbation of Kf^^, i.e., the 
limiting Kll^iti). The rank of the perturbation depends 
on the number of disconnected components comprising 
the highly-diffusive region. This special limiting form 
of S{m) allows us to build a robust approximation of 
S{m)^^ by merely using solvers for by the help 

of the Sherman-Morrison- Woodbury formula. In Sec- 
tion [s] by using SPA, the asymptotic behaviour of sub- 
matrices is provided and the final convergence proof of 
Backs is based on that result. 

The remainder of the article is structured as fol- 
lows. In Section [6] we reveal the qualitative nature 
of the solution of the high-contrast diffusion equation 
and the resulting decoupling in the solution. This valu- 
able insight is provided by SPA. In Section [7] we high- 
light the implementation aspects of Backs- In addi- 
tion, we present how subdomain deflation gives rise to 
a dramatic computational saving due to reduction to 
a block diagonal system. Finally, in Section [8] we nu- 
merically demonstrate the simultaneous robustness of 
Backs with respect to magnitude of the coefficient 
contrast and the mesh size. We also provide compar- 
isons to CCMG preconditioner with different prolonga- 
tion operators and smoothers. 



2 The underlying PDE and the linear system 



We study the following high-contrast diffusion equa- 



tion. 



-V • (aVu) = / in 12 
u = on df2 



(1) 



where J7 C M'*, d — 2,3. The coefficient a{x) may 
vary over many orders of magnitude in an unstruc- 
tured way on f2. Problems with high-contrast coeffi- 
cients are ubiquitous in porous media ftow applications. 
Many examples of this kind arise in groundwater flow 
and oil reservoir simulation; see for example the com- 
prehensive overviews [T|ITT |[20p3] . Consequently, devel- 
opment of efficient solvers for high-contrast heteroge- 
neous media has been an active area of research, specif- 
ically in the setting of multiscale solvers [51 IT^[TH[T51 
[21] . In addition, the fictitious domain method and com- 
posite materials are also sources of rough coefficients; 
see the references in |16j . Important current applica- 
tions deal with composite materials whose components 
have nearly constant diffusivity, but vary by several or- 
ders of magnitude. In composite material applications, 
it is quite common to idealize the diffusivity by a piece- 
wise constant function. Likewise, we restrict the diffu- 
sion process to a binary regime (see Figure [l| in which 
the coefficient a is a piecewise constant function with 
the following values: 



a{x) 



TO 3> 1, X e Hh, 

1, X £ n^. 



Due to the atomistic structure of matter, the physical 
treatment of diffusion involves regular (C°°-) diffusiv- 
ity. Aksoylu and Beyer 4 showed that the idealization 
of diffusivity by piecewise constant coefficients is mean- 
ingful by showing a continuous dependence of the solu- 
tions on the diffusivity. 

Let us denote the linear system arising from the 
finite volume discretization by: 

K{m) x{m) = h. (2) 

Let Q be decomposed with respect to diffusivity value 
as 




Fig. 1 i7 = Qu U where Qu and Qj^ are high and low 
diffusivity regions, respectively. 



(3) 



where Qh and fi^ denote the high and low diffusivity 
regions, respectively. When m-dependence is explicitly 
stated and the discretization system ^ is decomposed 
with respect to (|3|, i.e., the magnitude of the coefficient 
values, we arrive at the following 2x2 block system: 



KHH{m) KHL{ni) 





xnim) 








XLirn) 







(4) 
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Note that all the subblocks in Q have m-dependence. 
This is the main difference from the FE discretization in 
which only the i7i?-subblock has m-dependence. Hence, 
the perturbation analysis of the FV discretization be- 
comes more involved than that of the FE discretization. 
The exact inverse of K can be written as: 



K-^{m) = 



Ihh 




- K 



fjij{m)KHL{m) 
III 



Kn^Hi^) 
S'-i(m) 

Ihh 
-KLH{rn)K]j]j{m) III 



(5) 



where Ihh and III denote the identity matrices of the 
appropriate dimension and the S{m) is the Schur com- 
plement of Kfjijim) in K{m) given by 

S{m) = KLLim) - KLH{m)K^'H{m)KHL{m). (6) 

Let Mhh denote the Neumann matrix correspond- 
ing to the pure Neumann problem for the Laplace op- 
erator on Qh- We write an important decomposition 
which will be used in the analysis to come: 



KHH{m) ^ mNnH + ^{m) 



(7) 



3 The cell-centered finite volume discretization 

We start by giving the cell-centered FV discretization 
for the 5-point stencil used in this article; for more de- 
tails see [1^ . Let T = Vij with i, j = 1, . . . , rt^/^ be the 
mesh and the control volume be defined by 

ViJ = [Xi-i/2,Xi+i/2] X [yj-l/2,2/j + l/2]- 

The FV scheme is constructed by integrating ([T|) over 
each control volume Vij, which yields 



Vj + l/2 

yj-1/2 
yj+i/2 

Vj-l/2 

Xi+l/2 

Xi-1/2 
3:i + l/2 



ai+i/2j Ux{xi+i/2,y) dy 
ai-i/2.,j Ux{xi^i/2,y) dy 
ai,i-i/2 Uy{x,yj_i/2)dx 

7 + 1/2 Uy{x,yj + i/2) dx 

f{x, y) dxdy. 
Thus, the discretization scheme is 

Fi+\I1,3 - Pi-l/2,] + Fi,j + l/2 - Fi,3-l/2 = hi,3fi,j, 



where h. 



over Vi.j, and where 



■I. J — kx y~ ky, and /^.j is the mean value of / 
ku 



ij + l/2 



kx 
kx 
ki, 



{u{x^+i,yj) - u{xi,yj)} , 
Tla - {u{x^,yj + l) - u{xi,yj)} , 



and 



2 a, 



a. 



(8) 



Note that and ha^ are the harmonic means of the 
diffusion coefficients for the adjacent control volumes 
in X and y directions respectively. The discretization 
formula can be written explicitly as follows: 



- ha,u{xi+i,yj) ~ha^_^u{xi_i,yj) 

+ (ha, +ha. + ~haj ^i)u{Xi, yj) 

- ha.u{x,,yj + i) -haj^M^i^yj-l) = ^ij/ij- 



(9) 



In our binary diffusivity regime, for notational con- 
venience, we denote the harmonic mean by 



2m 



1 



(10) 



The harmonic mean is used to ensure the continuity 
of the flux across the control volume interfaces. As a 
result, this flavor of finite volume discretization enjoys 
the desirable property of mass conservation. 

One can write the discretization matrix entries a 
priori due to the formula (|9|. Hence, in 2D, we explic- 
itly state each contribution of the off-diagonals to the 
diagonal entry values in the following: 



[K{m)]pp - (11) 

' m + m + m + m, p e Iqi , 

m + m + m + h„i, p G Jf2i and non-corner, 

m + m + km + hrm p G /r2i and corner, 

< 1 +1 +1 +hra, pe Fn^, 

1+1+1 +1, 1^2 and strictly interior, 

5, p G 7/72 and non-corner bdry, 

^6, p G lor^ and corner bdry. 



3.1 Comparison of flnite element and flnite volume 
discretizations on a ID example 

We explicitly provide the discretization matrix utilizing 
the FV method in ^ corresponding to ([ij . The domain 
is chosen to be J7 := (0,1) with the highly-diffusive 
island fin '■= (2/7,5/7). The mesh consists of 7 cells. 
The cells and cell-centers are denoted by , ■ • ■ , ^7 and 
xi, . . . , XT, respectively. See Figure [2] 
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I -XQ - I - X4- I -X2 - I - 2^1 - I - Xz~ I -X^ - I - X-J- I 



Fig. 2 (Top) The finite volume mesh where the cell-centers are 
denoted hy Xi, i = 1, . . . ,7 . (Bottom) The finite element mesh 
whore the vertex locations are denoted by x^, i = 1, . . . , 6. 



The corresponding submatrices in ^ are given be- 
low: 

KHH{m) 
KHL{m) 

KLL{m) 



2m 


-777 


— 777 




— 771 777 









—m 


771 + hm 







0" 











- 


h-,n 0_ 




l+~hm 





-1 ' 





1 + h„i 


-1 


-1 





3 





-1 


3 



Moreover, from (j7|, we obtain 



A{m) = 



2777 —777 —777 

— 777 777 

— 771 777 



7^™_0 

Kr 



We readily see that the 777-dependence of the matrices 
KLH{rn), KHLi'iTT'), Kll^tti), and Z\(777) is eliminated 
asymptotically. For instance, 

Z\(777) = Z\°° + 0(777-l), 

"0 0] 

2 +0(777-!). 
2 



By maintaining a similar configuration (see Figure 
[2]) used for the FV case, the corresponding submatrices 
in ^ for the linear FE discretization are given below. 
DOF that lie on the interface are always included to the 
DOF of the highly-diffusive island. Hence, only Kf^j^ 
subblock has r77-dependence. 



2m 

—777 
—777 







— 777 — 777 
2777 —777 

777 + 1 

— 777 777 +1 



_ 

2 0' 
2 



A 



FE 



2m —777 —777 
—777 2771 —777 
—777 777 
—777 777 

0' 



10 

1 



4 Spectral analysis of the diagonally scaled 
finite volume discretization matrix 

Let {1, . . . , s} denote the index of islands in the domain. 
Let Nk, k — 1, . . . , s be the FV discretization matrix 
of — zi on the k-th island f2k with respect to 7^ with 
homogeneous pure Neumann boundary condition. Let 
C denote the set of all DOF in J7. We will analyze the 
behavior of the spectrum of the symmetrically scaled 
discretization matrix 

^(777) (diagi<r(777))~^/^ i^(777) (diagX(777))-!^^ . (12) 

We start by examining the dependence of the (p, q)-th 
entry of K(m) on 777, for each p S C. Let i^s+i denote 
the outer region of islands, i.e.. 



s+l 



k=l. 



We denote the cell-centers that are adjacent to fis+i 
and the ones that are interior to i?^, fc = 1, . . . , s by 
Fq^ and Iq^ , respectively. On the other hand, the cell- 
centers in the outer region Qg+i that are adjacent to 
the islands Qk, k = 1, . . . , s are denoted by /^f2a+i . 

We define the following index set for p,q (z C with 

1 k, p and q G ^2^., fc = 1, . . . , s 
[s + l, p or q£ Hs+i, 



Ap a q) 



k, 



l,fc, 
1, 



P^ I{2k, k'- 
P^Tn^, k 
p e f2s+i- 



1. 
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Also, we define X{p /\p) — T{p)- 

Note that the discretization matrix and its entries 
can be written as follows: 

s 

K{m) = J2mNk{l) + Ns+i{m), 
fc=i 

[K{m)]pq= ^ ae{m)[Ni]pq, 
where 



m, £ = 1, . . . , s 

1, e=s + i, 



(13) 



and by abuse of notation, we have defined Ng :— A^^ (1) 
for f = 1, . . . , s and N^+i := Ns+i{m). Then, 

[A{m)]pg 

1-1/2 
ai{m)[Ni]pq 
eex{pAq) 

-1/2 

jei{q) 

For p £ C, [A{m)]pq — ii p,q are not adjacent cell- 
centers and [A(m)]pp = 1. Furthermore note that, [A(m)]pq 
is m-dependent only if either p or q E Ua;=i s+i ■ 
It is sufficient to study only the single island case 
because single island expression for 

K{7n) = mNi + N2{m) (14) 

can be simply generalized to 

s 

K{m) ^miVfc + A,+i(m). 

k=l 



4.1 Perturbation expansion analysis for the upper 
bound of the smallest eigenvalue 

We devise a perturbation expansion analysis based on 
m, in order to study the m-dependent spectral behavior 
oi A(m). Hence, only the matrix entries A^rnjpq, p ^ q, 
that have m-dependence are considered where cells p 



and q are adjacent. Combining (111 and (14 1, one can 
deduce that 



(15) 



We will use ([15| in the below analysis. For clarity, we 
treat the perturbation expansion in full detail for the 
first case. 



Case 1: p E J^?! and q E f22, 
[A{m)]pq 

- {m[Ni]pp + [N2im)]pp}-^^^ [A^2(™)]p, {[iV2(m)], 
= {[N,]pp + m-i[^2(m)]pp}"'/' [^^2(^)1^, 

X {[N2{m)]qq}-'/' 
= m-'/' { ^-'m^p'/' ~ l/2[iVi];p3/2 [N2{fn)]pp 

+ 0{m-^)}{mm)]pq}{[N2{m)]qq}-'/' 

x{m-'[N^]qq + 0{vi-')y'^' 

Case 2: pe T^^^ and q E Fq^ , 
[A{m)]pq 

- {m[Ni]pp + [iV2(TO)]pp}"'/' [mNi]pq 

X {m[N,]qq + [N2{m)]qqr'/^ 

= {[Ni]pp + m-i[iV2(™)]pp}"'^' [Ni]pq 

X {[m]qq + m-'[N2{m)]qqy'^^ 

= m;p'/'mpq[N,]-q'/' + oim-') 

Case 3: p E Fn^ and q E loi , 

[A{m)]pq 

= {m[N,]pp + [N2{m)]pp}-'/^ [mN,]pq {m[N,]qq}-^^^ 

- {[Ni]pp + m-^[N2{m)]pp}-'^' [N,]pq {[N,]qq}-'/^ 
^ m-^'^'[N,]pq[N,]-^'/^ + 0{m'') 

Case 4: p E Fq^ and q E (^2, 
[A{m)]pq 

= {[N2{m)]pp}-'/' [N2{m)]pq {[N2{m)]qq}-'/' 
= {m-^[N2{m)]ppy^^^ m-^[N2{m)]pq 

X {m-i[A^2(m)],,}"'/' 
= {m-^[N^%p + 0(m-2)}-i/2 {^-i[^oo]^^ ^ Q^^- 

x{m-'[N^]qq + 0{m-')y'^' 

Remark 1 Numerically we observe that the smallest eigen 
value is 0{m^^); see Table[l] In the above analysis, all 
the cases except Case 1 give the same 0{m^^). Case 1 
result of 0(771"^/^) estimate is an artifact of the per- 
turbation expansion. The same artifact has appeared in 
the FE analysis; see |T3l equ. (5.14)]. 

We will use following modification of N2 for our fur- 
ther analysis: 

'0 iipEFn, 



[N^]pq otherwise. 



(16) 
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Table 1 A 2D example showing the condition numbers and eigenvalues of the finite volume discretization matrix K{m) and its 
diagonally scaled version Aim) in which the eigenvalues are sorted in ascending order. 



m 


K(K{m)) 


Xi{K{ra)) 




A6i(X(m)) 




\&4{K{m)) 


re(A(m)) 


Ai(A(m)) 




M{A{m)) 




A64(A(m)) 


10° 


2.627 X 10^ 


3.045 X 10" 


1 


7.695 X 10" 


7.848 X 10° 


8.000 X 10° 


2.553 X 10^ 


7.538 X 10" 


-2 


1.800 X 10" 


1 


1.925 X 10° 


102 


1.248 X 10^ 


3.237 X 10- 


1 


7.995 X 10" 


2.040 X 10^ 


4.040 X 10^ 


3.448 X 10^ 


5.784 X 10" 


-3 


1.362 X 10" 


1 


1.994 X 10° 


lO* 


1.235 X 10^ 


3.240 X 10" 


1 


8.000 X 10° 


2.000 X 10* 


4.000 X 10"* 


3.258 X 10"' 


6.139 X 10" 


-5 


1.346 X 10" 


1 


1.999 X 10° 


10^ 


1.235 X 10'' 


3.240 X 10" 


1 


8.000 X 10° 


2.000 X 10^* 


4.000 X 10^ 


3.256 X 10** 


6.142 X 10" 


-7 


1.346 X 10" 


1 


2.000 X 10° 


108 


1.235 X 10^ 


3.240 X 10" 


1 


8.000 X 10° 


2.000 X 10** 


4.000 X 10® 


3.256 X 10* 


6.143 X 10" 


-9 


1.346 X 10" 


1 


2.000 X 10° 




1.235 X 10" 


3.240 X 10" 


1 


8.000 X 10° 


2.000 X 10^° 


4.000 X lO^" 


3.256 X 10^° 


6.143 X 10" 


-11 


1.346 X 10" 


1 


2.000 X 10° 



Consider the reduced version of ( 14 1 
k{7n) = mNi + N2: 

and let A{m) denote the diagonally scaled version of 
K^m). Then m-independent A(m) has a single zero 
eigenvalue. Next, we proceed with the elementwise anal- 
ysis of A{m) — A{m). 
Case 1: p G F^^ and q G 



[A{m)], 



-1/2 



-1/2 



-1/2 



2 \qq 



-fO(m'3/2) 



and from (16 1 we get 



[A{m%, - 0. 
Therefore, 

[A{m%^ ~ [i(m)]p, = 0(m-i/2) 
Case 2: j3 G F^^ and q E f2i, 



and from (16 1 we get 



Thus, 

[A{m)]pg - [i(m)]p, = O(m-i) 
Case 3: p E Fq^ {nodes of f22), 
[A{m)]pq 



-1/2 



2 Ipq [^2°]qq^'^ ' 



and from (16 1 we get 



Thus, 

[Aim)]pg ~ [A{m)]pq - 0{m-^). 
Case 4: Otherwise, 
[A{m)]pg - [i(m)]p, = 0. 

Therefore, Amax(^(m) - A{rn)) = Oim'^/"^). 



Lemma 1 Let G and H be symmetric matrices of di- 
mension n X n. Then, for k — 1, . . . ,n, the following 
holds: 

Afe(G) + < Afc(G + H)< Afc(G) + A„,ax(-ff). 

Proof The result follows from Courant-Fischer mini- 
max Theorem; see [12j Corollary 8.1.3]. □ 

From Lemma [l] we have 

Amin 

{A{m)) < A,ni„(A(m)) + Amax(^(m) - A{m)) 

= A,nax(^(TO) - Mm)) 

Moreover, we have for all fc > 1, 

Xk{A{m)) > Afe(i(m)) -f A„ii„(yl(m) - A{m)) 
>Afc(i(l))-0(m-i)>C 

for some constant C independent of to, asymptotically. 
Thus, A(m) has a single eigenvalue approaching to zero 
while the rest eigenvalues is bounded away from 0. 



4.2 Lower bound for the smallest eigenvalue 

We aim to show the following lower bound for the small- 
est eigenvalue: 



X^in{A{m)) >C m ^. 



(17) 



For that, we will establish the below main estimates in 
the discussion to follow: 



K{m)x > K{l)x > m~^x^ K{m)x 



(18) 



Remark 2 Establishing the estimate ([18|) is not as sim- 
ple as in the FE discretization case due to TO-dependence 



of N2 in ( 14 1 . This requires further detailed matrix anal- 
ysis. 
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4-2.1 K{m)x > x'^K{l)x estimate 



For the decomposition in (14), it is straightforward to 



see that mNi > Ni for m > 1. Hence, in order to 



establish ( 18 1 , we concentrate on the following auxiliary 
estimate: 



x'^N2{m)x > x'^N2{l)x. 



(19) 



First, note that N2{m) has positive diagonal and 
negative off-diagonal entries. In addition, due to the 
discretization formula ([9|, it has a row sum equal to 
zero with the exception that the row sums correspond- 
ing to cell-centers that are adjacent to the boundary 
are positive. Hence, N2{m) is a diagonally dominant 
matrix. We further decompose N2{m) as follows: 



N2{m) = N2{m) + R2 



(20) 



where N2{m) is a symmetric matrix with positive di- 
agonal entries and a row sum equal to zero, and R2 := 
N2{m) — N2{m) is the remainder matrix. 

Lemma 2 Let G be a symmetric matrix and have a 
row sum equal to zero. In addition, let [G]pp > 0, then 
G is symmetric positive semi-definite (SPSD). 

Proof Let Ap be the p-th eigenvalue of G. Using Ger- 
schgorin's Theorem with the fact that G has a row sum 
equal to zero yields: 



< 



The result follows from < Ap < 2[G]pp. □ 

By Lemma [2] N2{ni) is immediately SPSD. In ad- 
dition, i?2 is a diagonally dominant matrix with non- 
negative diagonal entries, again by Lemma |2] R2 is 
SPSD. Now, we can conclude that N2{m) is SPSD. 



From (111, one observes that [K{m)]pp is monoton- 



ically increasing in m. This important property implies 
that 



[N2{m)]pp > [N2{1)] 



m > 1. 



(21) 



We need the following additional decomposition: 
N2{m)^N2{l) + R2{m). (22) 



Combining ([2T| and (22 1, we obtain [R2{'m)\pp > 0. 



Hence, Lemma [2] implies that R2{rn) is SPSD. Using 
this, we now arrive at the auxiliary estimate ( 19 1 in the 
following: 

x'^N2{m)x = x'^N2{l)x + x'^R2{m)x + x'^ R2X 
> x'^ N2{l)x + x'^ R2X 

= X^N2{1)X. 



Consequently, 

mx'^ K{m)x = mx'^ Nix + x'^ N2{ra)x 

> x'^NiX + x'^N2{m)x 

> x'^ Nix + x"^ N2{l)x 
= x^K{l)x. 

4.2.2 mx'^ K{l)x > x'^K{m)x estimate 



Combining ( 14 1 and ( 20 1 , we obtain 



K{m) = mNi + N2{m) + R2, 

where Ni and R2 are SPSD and independent of m, 
which yields the following for m > 1: 

mx^ Nix > x'^Nix, 
mx'^ R2X > x"^ R2X, 

Thus, it is sufficient to establish the auxiliary estimate 



mx N2{l)x > x'^N2{m)x, 



(23) 



in order to establish (18 1 



By using Remark [TT] one can also show that 
m[N2{l)]pp > [^2(m)]pp. (24) 
We will use the following decomposition: 
miV2(l) = N2{m) + R2{m), (25) 
where i?2(w) is the SPSD remainder matrix. Combin- 



ing (24 1 and (25 1, we obtain [i?2(»Ti)]pp > 0, which leads 



us to the auxiliary estimate (23 1. 
Hence, 

mx'^ K{l)x — mx^ Nix + mx'^ N2{l)x 
> mx'^Nix + x'^ N2{m)x 
= x'^ K{m)x. 

In conclusion, we have obtained the two main es- 



timates in (18). These yield the similar estimates for 
diag K{m). 

x'^diag K{m)x > x"^diag K{l)x 

> TO^^x-^diag K{m)x (26) 



From ([T8| and (|26|, and letting Ci := Aniin(^(l)), we 
get: 

x'^K{m)x > x'^K{l)x > Cix^diag K{l)x 
> Ci m^^x^diag K{m)x, 

yielding 

Ci < A,ni„(A(m)) < C2 mT^^^ ■ 
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Remark 3 The merit of diagonal scaling becomes ap- 
parent after studying the spectral behavior of A(m). 
We observe that the number of the eigenvalues of A 
of 0{m^^) depends only on the number of isolated is- 
lands. On the other hand, the number of the eigenval- 
ues of K{m) of 0{m) depends on the number of DOF 
of the islands, asymptotically. In the example we report 
in Table [l] K{m) has 3 eigenvalues approaching to oo 
and 61 bounded eigenvalues, whereas, A{'m) has only 
one eigenvalue approaching to zero. The reduction in 
the number of m-dependent eigenvalues is a desirable 
feature for fast convergence of Krylov subspace solvers. 



5 Singular perturbation analysis on matrix 
entries 

5.1 Preliminaries on matrix properties 

The discretization formula ([9| with the harmonic means 
(|8| gives rise to the matrix entries given in (11). Due to 



harmonic means, the m-dependence of the matrix en- 
tries corresponding to HL and LL couplings is asymp- 
totically eliminated. As a result, the submatrices K}ji^(m) 
and KLL^m) do not have m-dependence asymptoti- 
cally: 



KLH{m)=K " 
KLL{rn) 



(27) 



The analysis in this article mainly relies on this crucial 
property. 

To analyze the to- robustness of preconditioners based 
on ([s]) , we need to analyze the asymptotic behaviour of 
the block components Knui'm)^^ , S{m)^^ , and 
KLH{rn)KHH{m)~^ as to — > cx3. This is the purpose 
of Lemma [3] below. To prepare for this, we further de- 
compose DOF associated with Qh into a set of interior 
DOF associated with index / and boundary DOF with 
index F. This leads to the following further block rep- 
resentation of 



KHH{m) 



Kii{m) Kir{m) 
Kri{m) Krr{m) 



(28) 



By using (111, we can write a more refined expression 



for the block Krrifn)'- 
Krr{m) = icj^'^V) + h^of)., 
with 



D^^l diag(0, 



,0,1,, 



1,2,...,,2), 



(29) 



where the number of 0, 1, and 2 entries is equal to the 
cardinality of interior, non-corner interface, and corner 



interface cell-centers, respectively. Thus, we can char- 
acterize the Neumann matrix Mhh on Qu as described 
in ^ as follows: 

KHH{m) = mAfHH + ^{m) , (30) 
1 

(31) 



A{m) 



o4^; 

+ 0(to-1), 



(32) 



where A{m) is a diagonal matrix due to (29 1. Nhh is 
a SPSD matrix with a simple zero eigenvalue and asso- 
ciated constant eigenvector. If nu denotes the number 
of DOF in Qhi a suitable normalized eigenvector is the 

1/2 

constant vector with entries , which we denote by 

ch- We further write in block form as = {ej , ej,). 

Remark 4 The Neumann matrix in the FV discretiza- 
tion is completely different than that of the FE case. 
However, as expected, the simple zero eigenvalue and 
the associated constant eigenvector are maintained be- 
cause this is an inherent property of the underlying 



PDE. In addition, Z\ in (30) in the FE discretization is 



not necessarily diagonal and does not have TO-dependence. 

Finally we note that the ofl-diagonal blocks in ^ 
have the decomposition: 



KLH{m) = [OKLrim)] - KHLimY' 



(33) 



As we prepare for the proof of our main Lemma, we 
need to define the following quantity which will also be 



used for subdomain defiation in ( 50 1 
r]{m) := e^KHH{m)eH- 



(34) 



r]{m) > because KHH{m) is SPD as a diagonal sub- 



can reduce the expression in ( 34 1 to 



block of if (to). Moreover, combining (30) and (31), one 

(35) 



In fact, (35) can be expressed explicitly by: 



ri{m) = hr, 
where 



(36) 



nn nc '-^ ^{non-corner interface cell-centers} 
"fT-H c '■= #{corner interface cell-centers}. 



Finally, (|36| delivers the asymptotic convergence of ex- 

(37) 



pression for rj^m): 
r]{m) — rj + 0{m^^), 
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5.2 The main results on the preconditioner 

Lemma 3 The asymptotic behaviour of the submatri- 
ces in (|5| is described by the following: 

(i) KHH{m)-^ = CHir^el + 0{m-^), 
(ii) S{m) = Kf^ - {Kfrer)r)-\elKf^) + 0{m-^), 
(Hi) KLH{m)KHH{m)-^ = {Kf^rer)v~'el + 0{m-'). 

Proof Since Mhh is symmetric positive semidefinite we 
have the eigenvalue decomposition: 



Z^NhhZ = diag(Ai, . . . , A„^_i,0), 



(38) 



where {A^ : z = 1, . . . , njj } is a non-increasing sequence 
of eigenvalues ol Nhh and Z is orthogonal. Since, the 
eigenvector corresponding to the zero eigenvalue is con- 
stant, we can write Z 



we have 



Z I en 



and so, using (30 1 



Z'^KHH{m)Z 

m diag(Ai, . . . , A„„_ J + Z^ A{m)Z Z'^A{m)eH 
elA{m)Z elA{m)eH 

yl(m) 6{m) 



S (m) rj{m) 



(39) 



To find the limiting form of Khh{itt-) ^ note that 

A{m) = m diag(Ai, . . . , A„„_ J + Z'^A{m)Z 
= m diag(Ai, . . . , A„^_J 

X (/ + m-i diag(Ar\...,A-i_JZ^Z\(m)Z 



Now, let Cx := maxi<„^ A^ ^ and Ca be the con- 



stant in (321, i.e., ||Z\(m)||2 < \\A°^\\2 + CAm-\ Then, 



for sufficiently large m, 
\\A{mr%< 



m-^Cx 



< 



l-m-^Cx\\ZTA{m)Z\\2 



(40) 



where 



Ca ■■=^-T^'^Cx\\z^\\2\\A°^h\\z\\2 

-m-'' CACx\\~Z^U\Zh 

= 1+0(771^1). 



Hence, combining ( |40[ ) and (41 1, we obtain: 



(41) 
(42) 



We proceed with the following inversion: 
1 



A{m) 6{m) 
6{m)^ r}{m) 



where 

U{m) := 

V{m) := 



i -A{m)'~^5{jn) 
0^ 1 

i(m)-i 



0"^ [rjim) — 5{rii)^ Aim) ^5{m)^ 



Using ( 32 1 , one gets 



~6{m) = Z^A°°eH + 0{m-^). 



(43) 



Then, (|37f , (|42|), and (|43f imply that 

r/(m)~i = T]-^ + 0{m-^), 
[/(to) = l + 0{m-^), 

O 

0^ 

Combining the above results, we arrive at 



V{m) 



+ C'(m~i). 



A{m) 5{m) 
5(rn)^ rj{rn) 



-, -1 







0^ ?7-i 



and, by ( p39| , we have 
KHH{m)-^ = Z 



O 

0^ 



+ 0(to-1) , (44) 



Z'^ + ©(m-i) 



which proves part (i) of the Lemma. 

Parts (ii) and (iii) follow from simple substitution, 



using (|6| and (33 1 



□ 



of the proposed preconditioner ( 47 ) : 



We use the following limiting forms in the definition 

(45) 
(46) 



1. Then, 


^-oot 


-1 T 




QTh '■— 


^LH^HH-> 




S°° := 





Based on the above perturbation analysis we pro- 
pose the following preconditioner: 

KHHim)'^ 
S°° 



B 



AGKS 



(m) 



Ihh —Q 
III 



OO 

LH 



(47) 



'Ihh ■ 

The following theorem shows that Backs is an ef- 
fective preconditioner for to 3> 1 . 



Theorem 1 For to sufficiently large we have 
a{BAGKs{m) K{m)) C [1 - cm-^/^ 1 + cto-^/^j 
for some constant c independent of to, and therefore 
^{BAGKs{ni) K{m)) - 1 + 0{m-^'''). 
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Fig. 3 (Left) The matrix in (|49| corresponds to a homogeneous 
Dirichlet for the Laplacian on Q under the constraint that the so- 
lution is constant on Qh- (Right) xh^ (mi) = ch^ +0{m~ ), i = 
1,2 where E^nd correspond to square and triangle 

shaped highly-diffusive islands, respectively with rrii = 10^. 



Proof The TO-dependence of Kh l (™) , Klh [ttl) , and 
KLhitn) are ehminated asymptotically as in ( p?] ). Fur- 
thermore, m-dependence of rjirn) is also eliminated as 



a result of (32). Therefore, the result follows by slightly 



modifying the proof, i.e., by incorporating the m- de- 
pendencies of the above entities, given in Aksoylu et 
al. [7] for the FE case. □ 



6 Qualitative nature of the solution of the 
high-contrast diffusion equation and decoupling 

We advocate the usage of SPA because it is a very effec- 
tive tool in gaining qualitative insight about the asymp- 
totic behavior of the solution of the underlying PDE. 
Through SPA, in Lemma[3] we were able to fully reveal 
the asymptotic behaviour of the submatrices of K in 
(|5|. This information leads to a characterization of the 
limit of the underlying discretized inverse operator and 
we studied this in more detail in [5]. We now prove that 
asymptotically, the solution over the highly- diffusive is- 
land goes to a constant vector. This is probably the 
most fundamental qualitative feature of the solution of 
the high-contrast diffusion equation. 



Lemma 4 Let en be the constant vector. Then, 



Proof We prove the result by providing an explicit quan- 
tification of the limiting process based on Lemma |3] 



XL{m) = S l(m) {bL - Kj^Hy,,,j^yfj^y,,o)i 

= S^'\bL-K^HKfHbH} + 0{fr 
= : xf +0{m-^), 
xnirn) = K^\{7n) {bn - KHL{rn)xL{m)} 
= K^'H{bH-K^Lxf} + 0{m-^) 
=: CH en + 0{m-^). 



□ 



SPA helps to reveal a further qualitative property, 
namely, the decoupling phenomenon. We show how the 
original algebraic solution strategy decouples into two 
algebraic problems of different nature and the decou- 



pling is indeed an implication of ( 48 1 . This observation 



can be a promising research avenue when it is carried 
over to different classes of PDEs and discretizations. 
We will pursue this in our future research. 

In order to show the decoupling, let us start by not- 
ing that S°° in Lemma |3] can also be interpreted as the 
Schur complement of Killer in the matrix 



c KLF&r Kll 



for any nonzero value of c. In particular, if we choose 



'■H ' 



then cer — Ir, the vector of all ones on F 



and, using also ( 30 ) , we have 



Kll 



KlhI 



H 



Kll 



This is the stiffness matrix for a pure Dirichlet prob- 
lem for the Laplacian on all of J7 with the additional 
constraint that the solution is constant on fiu- See Fig- 
ure [3j Thus, when to 3> 1, the original problem decou- 
ples almost entirely into a (regularized) Neumann prob- 
lem (i.e. KfjH{m)) for the Laplacian on Qfj (scaled by 
m) and a Dirichlet problem (i.e. KK^^) for the Lapla- 
cian on all of J7, but under the additional constraint 
that the solution is constant on Qh- 

Next, we show an additional decoupling. This time 
the preconditioner decouples into a block diagonal ma- 
trix with the help of a deflation method. 



XH{m) = CH en + 0{m 



(48) 

7 Implementation aspects and the related 
deflation method 



where ch is determined by the solution in the lowly- 
diffusive region. 



The fact that Mhh has a simple zero eigenvalue (with 
the corresponding constant eigenvector en) and Mhh is 
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of co-rank 1 imply that KHni'm) has a single eigenvalue 
of 0(1) and n// — 1 eigenvalues of 0{m). For sufficiently 
large m, Ch can be taken as an approximate eigenvec- 
tor corresponding to JVhh'^ single smallest eigenvalue. 
Therefore, in the light of ([s]), 



[4,0^] 



becomes an approximate eigenvector corresponding to 
the smallest eigenvalue of the decoupled matrix: 



KHH{m) 
S{m) 



(49) 



In order to eliminate the negative effect of the small- 
est eigenvalue, we utilize a deflation method under the 



constraint that it gives the decoupling as in ( 49 1 . If such 



decoupling occurs, it would provide a large computa- 
tional advantage. We utilize a deflation method, known 
as suhdomain deflation P^[^E71E5] constructed by the 
following ii'-orthogonal projection onto the subspace 
orthogonal to e, provides the desired decoupling. 



/ — e ri{ra) ^ K . 



(50) 



We apply Backs within a conjugate gradient algo- 
rithm for the deflated system 



VKx^ = Vb, 



where x'^ := V^x is the projected solution. The com- 
ponent of X in the direction of e is then simply given 
by 



= {I — V'^)x= er]{m) 



-1 „T i 



When we rewrite ( 47 1 as 
BAGKsim) = D{m) L, 
where 



L 



Ihh 






5.00- 



D{m) := 

By observing the following interesting property 
LV^V, (52) 



the system in (51 1 after preconditioned by BACKsi'm) 
turns into: 

BAGKsim) (VKx^) = BAGKsim) (Vb) . 

Then, it reduces to the following block diagonal system: 

D{m) VKx^ ^ D{m) Vb. 



Remark 5 The fact that e becomes an approximate eigen- 
vector corresponding to the smallest eigenvalue of the 



decoupled matrix in ( 49 1 is not necessarily true for the 



original matrix K{m) because it is not block diagonal. 
Therefore, utilizing the above deflation method in the 
direction of e will not necessarily provide any further 
robustness for the underlying preconditioner if that pre- 
conditioner uses off-diagonal blocks of K{m) as well. 
Multigrid method is one such preconditioner. In fact, 
numerically we observed that introducing deflation did 
not improve the multigrid convergence rate at all. If one 
still wants to improve the convergence rate by the help 
of deflation, then the approximate eigenvector corre- 
sponding to the smallest eigenvalue must be computed 
in an alternative way rather than the simple usage of 
e. Consequently, we can say that Bagks{''tT')i by its 
design, naturally goes with subdomain deflation. The 
incorporation of subdomain deflation not only brings 
robustness with respect to the smallest eigenvalue, but 
also provides huge computational savings due to reduc- 
tion to a block diagonal system. This is a significant 
computational advantage Sacks (™) offers. 



By exploiting the fact that S°° in (46) is only a 
low-rank perturbation of K°° , we can build robust pre- 
conditioners for in ( 47 ) via standard multigrid pre- 



(51) conditioners. Combining (451 and (46 1, we arrive at 



VT] 



where v := K^j^en- If Mll denotes a standard multi- 
grid V-cycle for Kll, we can construct an efficient and 
robust preconditioner for S°° using the Sherman- 
Morrison- Woodbury formula, i.e. 



s- 



M 



LL 



LL- 



(53) 



Note also that we can precompute and store M^^v dur- 
ing the setup phase. This means we only need to apply 
the multigrid V-cycle Mll once per iteration. There- 
fore, the following practical version of preconditioner 



(47 1 is used in the implementation: 



Ba^ 



GKS 



Ihh -Q 
III 



00 

LH 



Mhh 

s-^ 



Ihh 
-Qfn III 



(54) 



8 Numerical experiments 

The goal of the numerical experiments is to compare 
the performance of the two preconditioners: AGKS and 
CCMG. The domain is unit square with a uniform mesh 
consisting of 2^ x 2^, ^ = 3, . . . ,6, cells. The coarsest 
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Table 2 Prcconditioner = CCMG, Prolongation = Wesseling-Khalil, Smoother = ILU 



h\m 


10" 


lOi 


102 


103 


10"' 10^ 10^5 10^ 10* lO'J 


1/8 


10,0.025 


29,0.472 


30,0.500 


35,0.528 


36,0.559 40,0.593 44,0.622 52,0.669 52,0.669 58,0.697 


1/16 


3,0.000 


3,0.000 


3,0.001 


6,0.024 


15,0.236 60+, 0.934 60+, 0.934 60+, 0.934 60+, 0.934 60+, 0.934 


1/32 


3,0.000 


3,0.000 


4,0.004 


9,0.074 


36,0.558 60+, 0.891 60+, 0.899 60+, 0.899 60+, 0.899 60+, 0.899 


1/64 


3,0.000 


4,0.001 


5,0.015 


9,0.084 


15,0.239 20,0.315 15,0.246 9,0.093 8,0.067 8,0.032 



Table 3 Preconditioner = CCMG, Prolongation = Wesseling-Khalil, Smoother = sGS 



h\m 


10° 


lOi 


10^ 


10^ 


10" 


10^ 10"^ 


10^ 


10« 


109 


1/8 


10,0.025 


29,0.472 


30,0.500 


35,0.528 


36,0.559 


40,0.593 44,0.622 


52,0.669 


52,0.669 


58,0.697 


1/16 


12,0.166 


14,0.222 


19,0.310 


26,0.409 


38,0.535 


52,0.649 60+, 0.777 60+, 0.843 60+, 0.938 


oo, 1.070 


1/32 


14,0.195 


15,0.217 


18,0.294 


25,0.432 


39,0.552 


58,0.698 60+, 0.917 


00,1.002 


oo, 1.080 


oo, 1.127 


1/64 


13,0.197 


14,0.221 


17,0.282 


22,0.362 


31,0.497 


48,0.645 60+, 0.793 


00,0.954 


oo, 1.097 


00,1.120 



level mesh contains 8x8 cells with a highly-diffusive 
single island of size 2x2 cells centered at the point 
(3/8,3/8). For the discussion of multiple disconnected 
islands, refer to [TJ Sections 3 and 4]. 

We denote the norm of the relative residual at iter- 
ation t by rr^*); 



r(*)| 



where r'^*^ denotes the residual at iteration t with a 
stopping criterion of rr^*) < 10~^. In the tables, we 
report the iteration count and the average reduction 
factor of the residual which is defined as: 



i/t 



We enforce an iteration bound of 60. If the the method 
seems to converge slightly beyond this bound, we de- 
note it by 60-I-. Whereas, the divergence is denoted by 
oo. In Tables [2j|9] iteration count and the average re- 
duction factor are reported for combinations of precon- 
ditioner, prolongation, and smoother types. 

We use Galerkin variational approach to construct 
the coarser level algebraic systems. There are two types 
of prolongation operators under consideration; Wesseling 
Khalil (3T] and bi-linear, given by respectively: 



piWK) ^ 



1 

16 



11 

13 2 

2 3 1 

11 

13 3 1 

3 9 9 3 

3 9 9 3 

13 3 1 



and RiWK)^piWKr 



2h 
h 



and R 



[B) ^ p(B) 



2h 



The multigrid preconditioner CCMG is derived from 
the implementation by Aksoylu, Bond, and Hoist [B]. 
We employ a V(l,l)-cycle with point symmetric Gauss- 
Seidel (sGS) or ILU smoothers. A direct solver is used 
for the coarsest level. We construct two different multi- 
level hierarchies for multigrid preconditioners Mjjh in 



(54 1 and M^^, in (53 1 for DOF corresponding to f^n 



and ]?!,, respectively. The corresponding prolongation 
matrices Phh and P^l are extracted from the prolon- 
gation matrix for whole domain Q in the fashion fol- 
lowing ([4|: 



P = 



Phh Phl 
Plh Pll 



The superior performance of AGKS preconditioner is 
partially due to employing these two distinct multi- 
level hierarchies, which is very much in spirit of the 
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Table 4 Prcconditioncr = CCMG, Prolongation = Bi-linear, Smoother = ILU 



h\m 


10° 


lOi 


10^ 


10^ 


10* 


10^ 


10"^ 


10^ 


10* 


10^ 


1/8 


10,0.025 


29,0.472 


30,0.500 


35,0.528 


36,0.559 


40,0.593 


44,0.622 


52,0.669 


52,0.669 


58,0.697 


1/16 


3,0.000 


3,0.000 


3,0.001 


6,0.003 


15,0.216 


00,0.954 


00,0.954 


00,0.954 


00,0.954 


00,0.954 


1/32 


3,0.000 


3,0.000 


4,0.003 


7,0.045 


19,0.316 


60+, 0.838 60+, 0.832 


60+, 0.842 60+, 0.843 


60+, 0.843 


1/64 


3,0.000 


3,0.001 


5,0.009 


7,0.028 


18,0.290 


16,0.239 


10,0.097 


8,0.050 


8,0.049 


8,0.056 



Table 5 Prcconditioncr = CCMG, Prolongation = Bi-lincar, Smoother = sGS 



h\m 


10° 


lOi 


102 


10^ 


10* 


10^ 10"^ 


10^ 


10* 


109 


1/8 


10,0.025 


29,0.472 


30,0.500 


35,0.528 


36,0.559 


40,0.593 44,0.622 


52,0.669 


52,0.669 


58,0.697 


1/16 


11,0.122 


13,0.194 


18,0.293 


24,0.414 


38,0.564 


58,0.692 60+, 0.827 60+, 0.931 


00,0.995 


oo, 1.094 


1/32 


8,0.075 


11,0.121 


15,0.229 


22,0.362 


35,0.545 


60+, 0.722 60+, 0.903 


00,1.022 


oo, 1.085 


oo, 1.133 


1/64 


8,0.068 


10,0.115 


14,0.216 


19,0.334 


27, 0.449 


44,0.600 60+, 0.773 


00,0.965 


oo, 1.116 


oo, 1.168 



aforementioned decoupling in Section [7] In fact, due to 
decoupling, AGKS technology allows the usage of any 
ordinary prolongation operator. This operator does not 
have to be constructed in a sophisticated manner as in 
the case of Wesseling and Khalil [31, or Kwak [17 . 

As emphasized in our preceding paper [7,, AGKS 
can be used purely as an algebraic preconditioner. There- 
fore, the standard multigrid preconditioner constraint 
that the coarsest level mesh resolves the boundary of 
the island is automatically eliminated. However, for a 
fair comparison, we enforce the coarsest level mesh to 
have that property. 

When the discretization matrix is scaled by 1/m, we 
observe a significant reduction in the iteration count for 
the AGKS preconditioner, while, CCMG suffers from 
inconsistent convergence behaviour. That is why, we 
only report the unsealed case for CCMG. Moreover, 
for the CCMG preconditioner, we use lexicographic or- 
dering. On the other hand, for AGKS, we follow the 
standard way of ordering the highly-diffusive after the 
lowly-diffusive DOF as used in [7]. 

Note that as the diffusion coefficient m increases, 
the CCMG method becomes less effective. For suffi- 
ciently large m, it even diverges. CCMG shows this 
adverse behaviour for all types of transfer operators, 
and smoother types for almost all levels. From Tables 
[2]-[5] CCMG preconditioner under ILU smoother pro- 
vides about three times faster convergence than that 
under sGS smoother. However, we also observe that the 



CCMG preconditioner becomes totally ineffective even 
diverges for m > 10^ under both types of smoothers; 
see the corresponding columns in Tables [2]-[5] 

On the other hand, the AGKS preconditioner be- 
comes more effective with increasing m. Even for the 
smallest m value, AGKS performance still remains com- 
parable to the CCMG performance. Furthermore, as 
seen in Tables |6}|9] AGKS preconditioner demonstrates 
consistently similar iteration counts and contraction num- 
bers for all types of transfer operators and smoothers. 
Therefore, AGKS performance is independent of the 
utilized prolongation operators and smoothers. 

CCMG performance heavily depends on the smoother 
choice. As opposed to expectations from the vertex- 
centered case, we observe that x- and y-line smoothers 
do not improve CCMG performance compared to sGS 
and ILU smoothers; see Llorente and Melson [12] for 
other ordering related smoother complications. As pointed 
out by Mohr and Wienands [5T], CCMG may need 
a sophisticated smoother like ILU. The specific ILU 
choice such as no-fill-in and sparsity threshold dramati- 
cally changes CCMG iteration counts. We use a sparsity 
threshold of 10"^ and observe that it improves the it- 
eration count by a factor of 3 compared to the no-fill-in 
case. The behaviour of ILU was extensively studied by 
Wittum [21133], also see [Ml pp. 98 and 134]. 

We observe an interesting cut-off m value for perfor- 
mances of preconditioners. While, CCMG performance 
starts to deteriorate at about m — 10^, the AGKS 
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Table 6 Preconditioner = AGKS, Prolongation = Wesseling-Khalil, Smoother = ILU 



h\m 


10° 


lOi 


10^ 10^ 


10-* 


10^ 


10^ 


lo'- 


10^ 


10^ 


10" 


10" 


1/8 


22,0.371 


10,0.115 


10,0.116 9,0.078 


9,0.056 


8,0.059 


8,0.045 


8,0.039 


8,0.039 


8,0.039 


8,0.039 


8,0.039 


1/16 


16,0.240 


13,0.199 


11,0.131 9,0.098 


8, 0.043 


7, 0.043 


6,0.018 


6,0.011 


6,0.007 


6,0.010 


5,0.010 


5,0.006 


1/32 


20,0.329 


19,0.313 


14,0.192 11,0.127 


9, 0.093 


8, 0.035 


7, 0.041 


6,0.018 


6,0.012 


6,0.008 


6, 0.004 


5,0.008 


1/64 


29,0.484 


26, 0.435 


17,0.284 12,0.161 10,0.092 


8,0.061 


8,0.026 


6,0.031 


6,0.016 


6,0.010 


6, 0.006 


5,0.013 



Table 7 Preconditioner = AGKS, Prolongation = Wesseling-Khalil, Smoother = sGS 



h\m 


10° 


IQi 


10^ 103 


10^ 


10^ 


10^ 


107 


108 


109 


10" 


10" 


1/8 


22,0.371 


10,0.115 


10,0.116 9,0.078 


9,0.056 


8,0.059 


8,0.045 


8,0.039 


8, 0.039 


8,0.039 


8,0.039 


8,0.039 


1/16 


16,0.268 


13,0.201 


11,0.131 9,0.098 


8, 0.043 


7, 0.043 


6,0.017 


6,0.010 


6,0.007 


6,0.004 


5,0.010 


5,0.005 


1/32 


20,0.350 


19,0.317 


14,0.192 11,0.127 


9, 0.093 


8, 0.035 


7, 0.041 


6,0.016 


6,0.010 


6,0.007 


6, 0.008 


5,0.008 


1/64 


29,0.483 


26, 0.436 


17,0.283 12,0.162 10,0.092 


8,0.061 


8, 0.030 


6,0.031 


6,0.017 


6,0.011 


6, 0.005 


5,0.013 



preconditioner reaches its peak performance and main- 
tains an optimal iteration count for all m > 10^. For 
instance, the iteration counts for the CCMG method 
jumps to 60 at about m = 10^ except for level four 
with ILU smoothing case. For that level, the iteration 
count starts to decrease although the numerical solution 
does not converge to the exact solution. Therefore, the 
performance of the CCMG preconditioner gets worse 
after m = 10^. Consequently, with respect to the mag- 
nitude of the coefficient contrast, this observation in- 
dicates that CCMG fails to be robust whereas AGKS 
maintains its robustness. 
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Table 8 Prcconditioncr = AGKS, Prolongation = Bi-linear, Smoother = ILU 



h\m 


10° 


lOi 


10^ 10^ 


lO* 


10^ 


10^ 


10^ 


10* 


10^ 


10" 


10" 


1/8 


22,0.371 


10,0.115 


10,0.116 9,0.078 


9,0.056 


8,0.059 


8,0.045 


8,0.039 


8,0.039 
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h\m 


10° 


lOi 


10^ 10^ 


10* 


10^ 


10^ 


10^ 


108 


109 


10" 


10" 


1/8 


22,0.371 


10,0.115 


10,0.116 9,0.078 


9,0.056 


8,0.059 


8,0.045 


8,0.039 


8,0.039 


8,0.039 


8,0.039 


8,0.039 


1/16 


16,0.256 


13,0.197 


11,0.131 9,0.098 


8,0.043 


7,0.043 


6,0.017 


6,0.010 


6,0.007 


6,0.004 


5,0.010 


5,0.004 
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20,0.352 
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